
close all;
clear all;
fclose all;

%create data structure
HCR=struct('CC',[],'tdTMLI',[],'tdTPC',[],'Globular',[]);


path1='C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\OTRcrex Ai14 Slc6a5 Aldh1a3 12.7.20\Processing\4. position measurements\';

dir_to_search=[path1];


% get the folder contents
d = dir(path1);
% remove all files (isdir property is 0)
dfolders = d([d(:).isdir]==1);
% remove '.' and '..'
dfolders = dfolders(~ismember({dfolders(:).name},{'.','..'}));

for SliceFolder=1:numel(dfolders) %slices
    temp_dir=[dfolders(SliceFolder).folder '\' dfolders(SliceFolder).name '\'];
    formatpattern = fullfile(temp_dir, '/*.mat');
    
    dInfo=dir(formatpattern);
    
    
    GolgiloadFileIndex=nan;
    CCloadFileIndex=nan;
    GlobularloadFileIndex=nan;
    MLI2loadFileIndex=nan;
    LugaroloadFileIndex=nan;
    
    for j=1:numel(dInfo) %Cell types
        %make logical for CC,globular etc in dInfo to load the correct file
        if ~isempty(strfind(dInfo(j).name,'CC'))
            CCloadFileIndex=j;
        else
            if   ~isempty(strfind(dInfo(j).name,'Globular'))
                GlobularloadFileIndex=j;
            else
                if ~isempty(strfind(dInfo(j).name,'tdTMLI'))
                    tdTMLIloadFileIndex=j;
                else
                    if ~isempty(strfind(dInfo(j).name,'tdTPC'))
                        tdTPCloadFileIndex=j;
                    else
                        
                    end
                end
            end
        end
    end
    
    %% Load CCs
    if ~isnan(CCloadFileIndex) %only run the section if there is a mat file for the corresponding Cell type, otherwise skip to next section
        
        load([dInfo(CCloadFileIndex).folder '\' dInfo(CCloadFileIndex).name]);
        HCR(SliceFolder).slice=CellsName(1:4);
        
        if exist('PC')
            HCR(SliceFolder).PC=PC;
            HCR(SliceFolder).ML=ML;
        end
        
        %calculate lobule length using PC layer boundaries and PC this is in
        %pixels of half resolution files of HCR from slidescanner
        if lobulePCboundaries(2)<lobulePCboundaries(1); % outline was traced in wrong direction, lobule X first.
            
            temp=[size(PC,1) lobulePCboundaries 1];
            temp=fliplr(temp); %need to be incremental for indexing purposes, will reverse again later after i'm done calculating length.
        else
            temp=[1 lobulePCboundaries size(PC,1)];
        end
        %
        
        for k=1:8
            lobuleSegments=PC([temp(k):temp(k+1)],:);
            d = diff(lobuleSegments);
            total_length(k) = sum(sqrt(sum(d.*d,2)));
        end
        
        if lobulePCboundaries(2)<lobulePCboundaries(1) % outline was traced in wrong direction, lobule X first reverse it again to get correct order
            total_length=fliplr(total_length);
        end
        
        HCR(SliceFolder).lobuleLength=total_length; % lobule I, II, IV-V,VI,VII,VII,IX,X
        
        
        
        
        % now onto CC stuff
        %asign lobule to specific cells
        lobuleID=[1 2 5 6 7 8 9 10];
        lobuleCellIdx=[1 lobuleCellIdx size(Cells,1)];
        Cellsperlobule=diff(lobuleCellIdx);
        HCR(SliceFolder).CCsperlobule=Cellsperlobule;
        
        clear lobule
        for l=1:8
            if ~(lobuleCellIdx(l)==lobuleCellIdx(l+1)) % add this in case there are no cells in some lobule (happens in OTRcre when looking at MLIs labeled with tdTom)
                
                lobule(lobuleCellIdx(l):lobuleCellIdx(l+1))=lobuleID(l)
            end
            
        end
        
        
        
        
        
        CCdata=[Cells distancePCL distanceML MLlength lobule']
        HCR(SliceFolder).CC=CCdata;
        clear Cells distancePCL distanceML MLlength lobule lobuleCellIdx CCdata
    end
    %% Load Globular
    if ~isnan(GlobularloadFileIndex) %only run the section if there is a mat file for the corresponding Cell type, otherwise skip to next section
        load([dInfo(GlobularloadFileIndex).folder '\' dInfo(GlobularloadFileIndex).name]);
        lobuleID=[1 2 5 6 7 8 9 10];
        lobuleCellIdx=[1 lobuleCellIdx size(Cells,1)]; %complete missing values (start and end).
        Cellsperlobule=diff(lobuleCellIdx); % calculate cells per lobule in order of lobuleID.
        HCR(SliceFolder).Globularperlobule=Cellsperlobule;
        clear lobule
        for l=1:8
            if ~(lobuleCellIdx(l)==lobuleCellIdx(l+1)) % add this in case there are no cells in some lobule (happens in OTRcre when looking at MLIs labeled with tdTom)
                
                lobule(lobuleCellIdx(l):lobuleCellIdx(l+1))=lobuleID(l)
            end
        end
        Globulardata=[Cells distancePCL distanceML MLlength lobule'];
        HCR(SliceFolder).Globular=Globulardata;
        clear Cells distancePCL distanceML MLlength lobule lobuleCellIdx Globulardata
    end
    %% load tdTPC
    if ~isnan(tdTPCloadFileIndex) %only run the section if there is a mat file for the corresponding Cell type, otherwise skip to next section
        load([dInfo(tdTPCloadFileIndex).folder '\' dInfo(tdTPCloadFileIndex).name]);
        lobuleID=[1 2 5 6 7 8 9 10];
        lobuleCellIdx=[1 lobuleCellIdx size(Cells,1)]; %complete missing values (start and end).
        Cellsperlobule=diff(lobuleCellIdx); % calculate cells per lobule in order of lobuleID.
        HCR(SliceFolder).tdTPCperlobule=Cellsperlobule;
        clear lobule
        for l=1:8
            if ~(lobuleCellIdx(l)==lobuleCellIdx(l+1)) % add this in case there are no cells in some lobule (happens in OTRcre when looking at MLIs labeled with tdTom)
                
                lobule(lobuleCellIdx(l):lobuleCellIdx(l+1))=lobuleID(l)
            end
        end
        tdTPCdata=[Cells distancePCL distanceML MLlength lobule'];
        HCR(SliceFolder).tdTPC=tdTPCdata;
        clear Cells distancePCL distanceML MLlength lobule lobuleCellIdx tdTPCdata Cellsperlobule
    end
    %% load tdTMLI
    if ~isnan(tdTMLIloadFileIndex) %only run the section if there is a mat file for the corresponding Cell type, otherwise skip to next section
        load([dInfo(tdTMLIloadFileIndex).folder '\' dInfo(tdTMLIloadFileIndex).name]);
        
        lobuleID=[1 2 5 6 7 8 9 10];
        lobuleCellIdx=[1 lobuleCellIdx size(Cells,1)]; %complete missing values (start and end).
        Cellsperlobule=diff(lobuleCellIdx); % calculate cells per lobule in order of lobuleID.
        HCR(SliceFolder).tdTMLIperlobule=Cellsperlobule;
        for l=1:8
            if ~(lobuleCellIdx(l)==lobuleCellIdx(l+1)) % add this in case there are no cells in some lobule (happens in OTRcre when looking at MLIs labeled with tdTom)
                
                lobule(lobuleCellIdx(l):lobuleCellIdx(l+1))=lobuleID(l)
            end
        end
        tdTMLIdata=[Cells distancePCL distanceML MLlength lobule'];
        HCR(SliceFolder).tdTMLI=tdTMLIdata;
        clear Cells distancePCL distanceML MLlength lobule lobuleCellIdx tdTMLIdata
    end
    
    
end

%% Plot shit
hold on
CCPCLdistance=[];
CCMLlength=[];
tdTMLIPCLdistance=[];
tdTMLIMLlength=[];
tdTPCPCLdistance=[];
tdTPCMLlength=[];
GolgiPCLdistance=[];
GolgiMLlength=[];
LugaroPCLdistance=[]
LugaroMLlength=[]

GlobularSliceCount=0;
tdTPCSliceCount=0;
tdTMLISliceCount=0;
CCSliceCount=0;
LugaroSliceCount=0;


for SliceFolder=1:numel(HCR)
    if  ~isempty(HCR(SliceFolder).CC)
        CCPCLdistance=[CCPCLdistance; HCR(SliceFolder).CC(:,3)];
        CCMLlength=[CCMLlength ; HCR(SliceFolder).CC(:,5)];
        CCSliceCount=CCSliceCount+1
    end
    
    if  ~isempty(HCR(SliceFolder).Globular)
        GlobularPCLdistance=[GlobularPCLdistance; HCR(SliceFolder).Globular(:,3)];
        GlobularMLlength=[GlobularMLlength ; HCR(SliceFolder).Globular(:,5)];
        GlobularSliceCount=GlobularSliceCount+1
    end
    
    if  ~isempty(HCR(SliceFolder).tdTMLI)
        tdTMLIPCLdistance=[tdTMLIPCLdistance; HCR(SliceFolder).tdTMLI(:,3)];
        tdTMLIMLlength=[tdTMLIMLlength ; HCR(SliceFolder).tdTMLI(:,5)];
        tdTMLISliceCount=tdTMLISliceCount+1
    end
    
    if  ~isempty(HCR(SliceFolder).tdTPC)
        tdTPCPCLdistance=[tdTPCPCLdistance; HCR(SliceFolder).tdTPC(:,3)];
        tdTPCMLlength=[tdTPCMLlength ; HCR(SliceFolder).tdTPC(:,5)];
        tdTPCSliceCount= tdTPCSliceCount+1
    end
    

end

binWidth=0.03
figure;CChist=histogram(CCPCLdistance./CCMLlength,'binwidth',binWidth,'FaceAlpha',0.5,'LineStyle', 'none')
figure;tdTMLIhist=histogram(tdTMLIPCLdistance./tdTMLIMLlength,'binwidth',binWidth,'FaceAlpha',0.5,'LineStyle', 'none')

% histogram(GlobularPCLdistance./GlobularMLlength,'binwidth',0.05,'DisplayStyle','stairs', 'EdgeColor', [0 0 0])
% histogram(GolgiPCLdistance./GolgiMLlength,'binwidth',0.05,'FaceAlpha',0.5,'LineStyle', 'none')
% histogram(MLI2PCLdistance./MLI2MLlength,'binwidth',0.05,'FaceAlpha',0.5,'LineStyle', 'none')



%%
color5=[27 192 218]./255;
color1=[0.29, 0.72, 0.35];
color2=[233 83 30]./255; % MLI2 color
color6=[255 0 35]./255; % OTRcre tdTomato cells color from illustrator
color7=[74 184 89]./255; %Globular color
color3=[0 0 0];
color4=[78 233 30]./255;



figure


CCStair=stairs(CChist.BinEdges(1:end-1), CChist.Values./CCSliceCount)
% Get the (y,x) coordinates from the original inputs or from
% the xdata and ydata properties of the stair object which will
% always be row vectors.
bottom = 0; %identify bottom; or use min(b.YData)
x = [CCStair.XData(1),repelem(CCStair.XData(2:end),2)];
y = [repelem(CCStair.YData(1:end-1),2),CCStair.YData(end)];
% plot(x,y,'y:') %should match stair plot
% Fill area
area(x,y,'FaceColor',color6,'Facealpha', 0.5,'EdgeColor',color6, 'EdgeAlpha',1,'LineWidth',2)



vline(0,'--k')
xlim([-1 1])
ylabel('count')

% MLI layer distribution

figure
tdTMLIStair=stairs(tdTMLIhist.BinEdges(1:end-1), tdTMLIhist.Values./tdTMLISliceCount)
% Get the (y,x) coordinates from the original inputs or from
% the xdata and ydata properties of the stair object which will
% always be row vectors.
bottom = 0; %identify bottom; or use min(b.YData)
x = [tdTMLIStair.XData(1),repelem(tdTMLIStair.XData(2:end),2)];
y = [repelem(tdTMLIStair.YData(1:end-1),2),tdTMLIStair.YData(end)];
% plot(x,y,'y:') %should match stair plot
% Fill area
area(x,y,'FaceColor',[255 165 0]./255,'Facealpha', 0.5,'EdgeColor',[255 165 0]./255, 'EdgeAlpha',1,'LineWidth',2)



vline(0,'--k')
xlim([-1 1])
ylabel('count')



















%% Make bar graph of density distribution

%make matrix for the numbers of each cell type (2d Vector: celltype(lobule, slice)
CC_lobules=nan(8,numel(HCR));


for SliceFolder=1:numel(HCR)
    CC_lobules(:,SliceFolder)=HCR(SliceFolder).CCsperlobule;
    tdTMLI_lobules(:,SliceFolder)=HCR(SliceFolder).tdTMLIperlobule;
    tdTPC_lobules(:,SliceFolder)=HCR(SliceFolder).tdTPCperlobule;
    
    
  
    lobule_length(:,SliceFolder)=HCR(SliceFolder).lobuleLength;
    
end

%% Do bar graphs first without any normalization
figure;
bar(nanmean(CC_lobules,2))
hold on;
for SliceFolder=1:numel(HCR)
    scatter([1:size(CC_lobules,1)] ,CC_lobules(:,SliceFolder),'filled')
   
end
box off
title('CCs')

figure;
bar(nanmean(tdTMLI_lobules,2))
hold on;
for SliceFolder=1:numel(HCR)
    scatter([1:size(tdTMLI_lobules,1)] ,tdTMLI_lobules(:,SliceFolder),'filled')
   
end
box off
title ('tdTMLI')

figure;
bar(nanmean(tdTPC_lobules,2))
hold on;
for SliceFolder=1:numel(HCR)
    scatter([1:size(tdTPC_lobules,1)] ,tdTPC_lobules(:,SliceFolder),'filled')
   
end
box off
title ('tdTPC')



%% make bar plots normalized to PC layer length
%analysis for length of PCL was done in Full resolution
%files for this set of HCR bc there was no registration. Full res scale is 160.67 nm/px
scale_factor=160.67/1e3; %convert it to microns/px.
lobule_lengthMicron=lobule_length*scale_factor

CC_lobulesNorm=CC_lobules./lobule_lengthMicron*100; %number of cells per 100 micron of PCL length
tdTMLI_lobulesNorm=tdTMLI_lobules./lobule_lengthMicron*100; %number of cells per 100 micron of PCL length
tdTPC_lobulesNorm=tdTPC_lobules./lobule_lengthMicron*100; %number of cells per 100 micron of PCL length


figure;
b=bar(nanmean(CC_lobulesNorm,2))
b.EdgeColor='none'
b.FaceColor=color6;
hold on;
for SliceFolder=1:numel(HCR)
    if isempty(strfind(HCR(SliceFolder).slice,'b3'))
        markerColor= [0.5 0.5 0.5];
    else
        markerColor=[0 0 0];
    end
    scatter([1:size(CC_lobulesNorm,1)] ,CC_lobulesNorm(:,SliceFolder),'filled','MarkerFaceColor', markerColor)
end
ylim([0 3])
title('CC norm')
box off
set(gcf,'color','none');
set(gca,'color','none');


figure;
b=bar(nanmean(tdTMLI_lobulesNorm,2))
b.EdgeColor='none'
b.FaceColor=[255 165 0]./255;
hold on;
for SliceFolder=1:numel(HCR)
    if isempty(strfind(HCR(SliceFolder).slice,'b3'))
        markerColor= [0.5 0.5 0.5];
    else
        markerColor=[0 0 0];
    end
    scatter([1:size(tdTMLI_lobulesNorm,1)] ,tdTMLI_lobulesNorm(:,SliceFolder),'filled','MarkerFaceColor', markerColor)
end
ylim([0 3])
title('tdTMLI norm')
box off
set(gcf,'color','none');
set(gca,'color','none');



figure;
b=bar(nanmean(tdTPC_lobulesNorm,2))
b.EdgeColor='none'
b.FaceColor=[210 105 30]./255;;
hold on;
for SliceFolder=1:numel(HCR)
    if isempty(strfind(HCR(SliceFolder).slice,'b3'))
        markerColor= [0.5 0.5 0.5];
    else
        markerColor=[0 0 0];
    end
    scatter([1:size(tdTPC_lobulesNorm,1)] ,tdTPC_lobulesNorm(:,SliceFolder),'filled','MarkerFaceColor', markerColor)
end
ylim([0 6])
title('tdTPC norm')
box off
set(gcf,'color','none');
set(gca,'color','none');



%% plot cell locations in all slices

CCscolor=color6;

alphaPCL=0.1;
alphaML=0.7;
mksize=2.5;
lobuleIdx=[1:8]
fig1=figure(1);
fig1.Renderer='Painters';
set(gcf,'position',[680 63 745 915])

for sliceNumber=1:numel(HCR)
    
    subaxis(numel(HCR),1,1+((sliceNumber-1)*1), 'SpacingVert',0,'spacingHoriz',0)
    %CCs
    Cells=HCR(sliceNumber).CC(:,[1:2]);
    PC=HCR(sliceNumber).PC;
    ML=HCR(sliceNumber).ML;
    
    % plot shit
    plot(PC(:,1),PC(:,2),'color', [0 0 0 alphaPCL]);
    hold on
    plot(ML(:,1),ML(:,2),'color', [0 0 0 alphaML]);
    temp=[1 cumsum(HCR(sliceNumber).CCsperlobule)+1];
    for l=2:numel(temp)
        
        scatter(Cells([temp(l-1):temp(l)],1),Cells([temp(l-1):temp(l)],2),mksize,'MarkerFaceColor',CCscolor,'MarkerEdgeColor','none','MarkerFaceAlpha',1);
    end
    xlim([0 25e3])
    box off
    set (gca,'Ydir','reverse')
    a=600
    ylim([0 25e3])
    axis square
    axis off
    
    
    
    
end


figure
tdTPCcolor=[210 105 30]./255;

alphaPCL=0.1;
alphaML=0.7;
mksize=2.5;
lobuleIdx=[1:8]
fig2=figure(2);
fig2.Renderer='Painters';
set(gcf,'position',[680 63 745 915])

for sliceNumber=1:numel(HCR)
    
    subaxis(numel(HCR),1,1+((sliceNumber-1)*1), 'SpacingVert',0,'spacingHoriz',0)
    %tdTPC
    Cells=HCR(sliceNumber).tdTPC(:,[1:2]);
    PC=HCR(sliceNumber).PC;
    ML=HCR(sliceNumber).ML;
    
    % plot shit
    plot(PC(:,1),PC(:,2),'color', [0 0 0 alphaPCL]);
    hold on
    plot(ML(:,1),ML(:,2),'color', [0 0 0 alphaML]);
    temp=[1 cumsum(HCR(sliceNumber).tdTPCperlobule)+1];
    for l=2:numel(temp)
        
        scatter(Cells([temp(l-1):temp(l)],1),Cells([temp(l-1):temp(l)],2),mksize,'MarkerFaceColor',tdTPCcolor,'MarkerEdgeColor','none','MarkerFaceAlpha',1);
    end
    xlim([0 25e3])
    box off
    set (gca,'Ydir','reverse')
    a=600
    ylim([0 25e3])
    axis square
    axis off
    
    
    
    
end



figure
tdTMLIcolor=[255 165 0]./255;

alphaPCL=0.1;
alphaML=0.7;
mksize=2.5;
lobuleIdx=[1:8]
fig3=figure(3);
fig3.Renderer='Painters';
set(gcf,'position',[680 63 745 915])

for sliceNumber=1:numel(HCR)
    
    subaxis(numel(HCR),1,1+((sliceNumber-1)*1), 'SpacingVert',0,'spacingHoriz',0)
    %tdTPC
    Cells=HCR(sliceNumber).tdTMLI(:,[1:2]);
    PC=HCR(sliceNumber).PC;
    ML=HCR(sliceNumber).ML;
    
    % plot shit
    plot(PC(:,1),PC(:,2),'color', [0 0 0 alphaPCL]);
    hold on
    plot(ML(:,1),ML(:,2),'color', [0 0 0 alphaML]);
    temp=[1 cumsum(HCR(sliceNumber).tdTMLIperlobule)+1];
    for l=2:numel(temp)
        
        scatter(Cells([temp(l-1):temp(l)],1),Cells([temp(l-1):temp(l)],2),mksize,'MarkerFaceColor',tdTMLIcolor,'MarkerEdgeColor','none','MarkerFaceAlpha',1);
    end
    xlim([0 25e3])
    box off
    set (gca,'Ydir','reverse')
    a=600
    ylim([0 25e3])
    axis square
    axis off
    
    
    
    
end





